function [R, llh] = expectation(X, model)
    mu = model.mu;
    Sigma = model.Sigma;
    w = model.w;                           
    
    n = size(X,2);                         
    k = size(mu,2);                        
    R = zeros(n,k);                        
    for i = 1:k                            
        R(:,i) = loggausspdf(X,mu(:,i),Sigma(:,:,i));
    end
    R = bsxfun(@plus,R,log(w));            
    T = logsumexp(R,2);                    
    llh = sum(T)/n; 
    R = exp(bsxfun(@minus,R,T));           
end

